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Eight isothermal equations of state are analyzed to yield quantitative measures of the degrees to 
whieh equation pairs ean be discriminated for real data, data of limited span and precision. Calculated 
curves allow one to assess the span and precision necessary in P-V data to allow unambiguous dis- 
crimination of various pairs. Some discussion is presented of bias and systematic error which may 
arise in least squares fitting. Usin<: exact synthetic data, we also illustrate for seven equation pairs the 
very large relative systematic errors in parameter and standard deviation estimates which arise 1 from 
such fitting of data of limited span with an incorrect but "close" equation model. General conclusions 
following from these results are discussed. Although the present work is principally concerned with 
discrimination between equations of state, its results are pertinent to the more general problem of 
choosing a "best" analytical model (linear or nonlinear) to represent experimental results. 
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1. Introduction 

Virtually all physical science is concerned at some 
stage with comparing experimental data with theo- 
retical predictions. Although no theories are ever fully 
verifiable, one nearly always wants to find that theo- 
retical model, from the limited set of possible models 
under consideration, which best represents the data, 
which allows the underlying phenomena to be better 
understood, and, if possible, which allows prediction 
outside the range of the original measurements. In the 
relatively early stages of investigation of a given 
domain, one usually does not know which of several 
theoretical or empirical models is likely to be most 
appropriate. This state of affairs is particularly likely 
to occur when the physical situation being studied is 
too complex to allow a tractable theoretical idealization, 
which is still sufficiently close to the experimental 
situation, to be accurate. Many-body interaction prob- 
lems, such as that of determining the exact equation of 
state of a solid or liquid, fall in this category. 

The problem of model discrimination is made dif- 
ficult by the presence of random and systematic errors 
in the data. In the present paper, it is assumed that 
systematic error in the data is absent or at least neg- 
ligible compared with other error. Systematic error can 
still be generated, of course, by the choice of an in- 
appropriate model [1|J and a question of considerable 
importance is: Under what conditions is it possible to 
discriminate adequately between several more-or-less- 
appropriate models, or equations? In the present paper, 
we shall be concerned with typical synthetic equation- 
of-state data generated without significant error of any 
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kind, reserving a detailed discussion of the effect of 
random errors to a later paper. It will be shown that 
by using such exact "data" we can investigate what 
sort of discrimination is possible between various 
equations of state in practical cases where measure- 
ments are of limited precision. 

In real life, experimental data have only limited 
accuracy and precision and always extend only over a 
limited range of the variables involved. This state of 
affairs suggests intuitively that one will be unable to 
discriminate adequately between two or more analyti- 
cal models which are sufficiently close together in their 
predictions for the range considered. We are here con- 
cerned with ways of making this intuition quantitative 
at least for the specific equations considered here. 
Since better discrimination may sometimes appear 
possible than is actually the case, just because of the 
presence of more or less random errors which happen 
to fall in a particular way, it is important to consider 
exact data before data with random errors. 

Although all that is often required of an equation of 
state, or more generally, a mathematical model of 
experimental results, is that it serve adequately as an 
interpolation and smoothing device for the data, the 
problem of model discrimination is usually still present 
even in this case. Unless the first model fitted passes 
all tests of adequacy, more than one model must be 
examined and a choice of available models made. The 
present paper discusses some general methods of 
model discrimination with specific illustrations taken 
from the equation of state field. Here we are concerned 
additionally with the task of estimating physically 
significant parameters of the material which led to the 
data in question. 

Two somewhat different situations frequently arise 
in the equation of state area. Often one starts with no, 



441 



or only crude, knowledge of the underlying parameters 
of the material under investigation. These parameters 
are then determined by fitting various equations of 
state to P-V data, usually by least squares techniques 
[1]. The most appropriate, or "best fit" equation will 
usually be that which leads to minimum estimated 
standard deviations of the fitted data points and of the 
parameters. The values of the parameters obtained 
from this fit are taken to be the best available estimates 
of the unknown material parameters. In general, how- 
ever, such values will not usually be good estimates 
unless the choice of model is appropriate for the data 
and leads to randomly distributed, essentially sto- 
chastically independent residuals, and the fitting 
procedure itself leads to negligibly biased parameter 
estimates. 

Sometimes one is able to obtain estimates of the 
parameters by other means than from fitting of direct 
P-V measurements. Now becoming popular for this 
purpose is the method of ultrasonic velocity measure- 
ments under pressure [2, 3|, pioneered by Lazarus 
[4J. Once having parameter estimates available for a 
certain material, one can, given the appropriate equa- 
tion of state, calculate volume values {or a range of 
pressure. Of course, with a limited number of param- 
eters available, as is always the case, calculated 
volumes can generally only be expected to remain 
reasonably accurate over a limited range of pressure. 
The nub of the problem here is usually in knowing 
what equation of state to use (and how far to trust it). 
Sometimes the determination of the best equation 
may be made by comparing ultrasonically derived 
parameters with those obtained from a least squares fit 
of direct P-V data for the same material. 

In either approach, one eventually obtains a set of 
parameters believed to be appropriate for the material 
under investigation. Although in actual practice these 
parameters will always be uncertain to some degree, 
it is nevertheless useful to ask, as a limiting case, how 
well one can distinguish between various equations of 
state when the parameters are actually exact (or are so 
considered) but when available P-V data are of limited 
precision. Some answers to this question are discussed 
later for eight different equations of state of some 
current interest. 

One of the important purposes of the present work is 
to point out that uniqueness is a limit seldom achieved 
in practice. Frequently an experimenter chooses a 
model to represent data of given range with the impli- 
cation or statement that the chosen model is "best" or 
"most applicable" without realizing or investigating 
sufficiently to find that other different models are 
equally applicable for the given data. 

Although the present analysis is concerned with 
discrimination between eight specific equations of 
state and thus involves quantitative results only for 
these equations, we expect that the results will also 
apply at least qualitatively to other not-too-different 
equations. More importantly, perhaps, the present 
discrimination methods and general approach can and 
should be applied to any experimental situation where 
it is important to establish one or more adequate 



mathematical representations of the data or, better, 
of the underlying process which led to the predictable 
part of the data. 

2. Equations of State Considered 



The material parameters with which we shall be 
concerned, all for isothermal conditions, are the 
specific volume, Vq, at a given reference pressure Po; 
the bulk modulus at P = Po, K {) = - Vo(dP/dV) \r=i> 
and various pressure derivatives of the bulk modulus, 
K, also evaluated at P=Po. For simplicity, let p = P— Po; 
then V=V {) at p = 0. Now Ki = r) = (dK/dP) | P =o, 
and K'd = (d 2 £/^P 2 ) | p=0 . The symbol rj has been 
introduced to simplify subsequent equations; it is 
dimensionless. It is also useful to introduce the further 
dimensionless quantity t// = KoK'd. Finally, define the 
dimensionless pressure variable z = p/Ko and the 
dimensionless density variable x = pjpo = Vo/V. 

Barsch and Chang [3] have recently given values for 
the parameters of Csl at 25 °C, plus temperature 
derivatives of these parameters. The quantity Vo may 
be calculated from x-ray measurements of the lattice 
constant. Other parameters such as /£<>, Kq, and K'i 
were obtained from ultrasonic measurements. Using 
the Barsch and Chang results, we have calculated the 
values of V {) , Ko, r\ and i// which then apply to Csl 
at 150 °C, an arbitrary choice of temperature. These 
values, as used in our computer studies, have 14 figure 
accuracy and may be considered the accurate values of 
some hypothetical material close to Csl at 150 °C. 
Of course as applied to Csl itself, only a few places 
in each parameter value are significant. To five 
figures, the parameter values are: V {) = 1.0184, 
K {) = 1.0503 X 10 2 kbar, n = 6.0382, and i// = -6.9897. 
Here we have taken Po = and Vo at 25 °C as unity. 
Thus, all volumes used here are reduced specific 
volumes and are dimensionless. The original Barsch 
and Chang 25 °C values are F =l, K = 118.9 ±0.6 
kbar, n = 5.86±0.11 and ^ = -0.052 ±0.002 kbar- 1 . 
These results lead to * = -6.2 at 25 °C. 

We shall be interested here in comparisons of, or 
discrimination between, eight different equations fre- 
quently employed in equation of state studies [1, 3J. 
We have adopted the approach of Barsch and Chang 
of designating the ordinary Murnaghan equation as the 
first-order Murnaghan equation (MEi), and the equa- 
tion previously [1J termed the second-order equation 
(SOE) as the second-order Murnaghan equation (ME 2). 
There are several forms of this latter equation, depend- 
ing upon the values of if- and ^P; here only one of these 
forms is pertinent. All eight of the equations are given 
in the form z=f(x) in table 1, which also lists acronyms 
for each equation. Some, but not all of them, may be 
expressed in inverse form, with x as an explicit function 
of z. Note that three of the equations are "first-order" 
in the sense of Barsch and Chang [3]. They involve 
n=S parameters: Vo, X () , and 77. The other five "second- 
order" equations involve ^ in addition. Finally, 
table 1 includes values of K ' = (dK/dP) _> x . 
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Table 1. Equations of state of interest written in the dimensionl ess form z=f(x) 



Equation 


Acronym 


Form 


n 


*: 


2 - P IK„=f(x) f(x) x - p/po - Vo/V 


Usual Tait 


UTE 


(ij+1) '[ ex|»{fr,+ l)(l-r')}-l] 


3 


— OC 


First-order 
Murna<jhan 


ME, 


J) '[.v"-l] 


3 


T? 


Second-order 
Murnaghan 


ME 2 

(•n 2 2*2t/0 


2{x (r <- -^ -l)/[(r/--2v|i) ,2 (x^'---^ ,2 + 1) 

-T)(X (1 >'~- 2 * )U2 -1) | 


4 


— 00 


Keane 


KE 

(—n 2 <i//<0) 


\T) :i l(-n 2 + ilj)-\[xW^>ln-\]-[iljl(r) 2 + 4j)\ \nx 


4 


Or+i/O/17 


First-order 
Birch 


BE, 


CM2)[x'"-x^][\ + (3/4)(t,-4)(*^- 1 ) ] 


3 


3 


Second-order 

Birch 


BE 2 


CV2)|x 7:; -x : ' :! J[M-(;V4)(r,-4)(x- ;{ -l) 
■f(l/24){143-r-9^(T)-7)4-9i/i}<a; 2 / :} -l) 2 ] 


4 


11/3 


Third-degree 


3SE 


(1-jr , ) + (1/2)(t;+1)(1-.v ' )- + (l/6)( if 


4 


-3 


Slater 




-f3T) + 2 + l/0 (1-jk ') :t 






Third-degree 
Davis-Gordon 


3DGE 


(v-l) + (1/21(7?- 1 )(.v- 1 ) 2 

+ (1/6)(t?--3t/ + 2 + «|;)(a:-1) ! 


4 


3 



The Keane equation only applies when — n 2 < ty < 0, 
conditions satisfied for the present parameter values. 
Although all of the equations must become poor 
models for sufficiently high z, failure is particularly 
evident for the UTE and MEj. The volume predicted 
by the UTE goes through zero at the finite z value of 
(77-f l) _1 [exp (17+ 1) — 1]. For the present form of the 
ME>, K' = U)K/c)P)=0 at z = -r)IV and V = at 
z = 2l[(r l ' 1 -2Vyi' 1 -r)]. The 3SE also suffers from the 
disadvantage that it predicts zero volume at finite 
pressure. All the other equations require infinite z to 
produce zero volume. 

The equations of table 1 are discussed in greater 
detail elsewhere [1, 3J. Although most of them have 
some macroscopic or phenomenological theoretical 
justification, here they may simply be regarded as 
empirical equations likely to be of some value in the 
P-Karea. 

3. Model Differences and AV Discrimination 

In order to examine differences in the predictions of 
the various models, we have, for a given set of /; or z 
values, calculated corresponding dimensionless V 
values, using in each equation the same 14-figure 
parameter values already mentioned. The V values 
were calculated using 14 figures, by iteration when 
necessary, with a resulting 13-figure or better accuracy. 
Finally, differences between V values of each possible 
pair of equations were calculated for each z value. The 
differences obtained for p = ll kbar, or 2 =0.1047, are 



listed in table 2, all multiplied by 10 4 for convenience. 
The A^s shown are formed by taking the V of one of 
the equations listed in the left column and subtracting 
from it the V calculated using one of the equations in 
the top row. Since the MEj— 3SE AV value is largest 
of all, the MEi yields the largest and the 3SE the 
smallest V value for this value of z. Similarly, we see 
that the BE2 and KE volume predictions are closest 
together here. 

In addition, in figures 1 to 5 we have plotted AV 
versus z for a variety of equation pairs. The boxed 
equation name is the equation from whose V values 
those of the equations named on the curves are sub- 
tracted. These five figures contain AV curves for most, 
but not quite all, of the possible pairs of equations. 
Curves have not been duplicated. Thus (F BE2 — ^bei) 
appears in figure 3 for BE 2 but not its negative in figure 
5 for BE 1. Negative values are indicated by using 
dashed lines. 

Table 2. Scaled volume differences, 10 4 AV . for ('(/no- 
tion pairs at z = 0.1047 



Equations 


3SE 


BE, 


ME, 


UTE 


3D(;e 


BE 2 


KE 


ME, 


3SE 



















BE, 


2.1 

















ME 2 


2.9 


0.77 















UTE 


3.7 


1.6 


0.81 













3DGE 


3.9 


1.8 


1.1 


0.27 











BE, 


4.2 


2.1 


1.3 


0. 19 


0.22 









KE 


4.3 


2.2 


1.4 


0.62 


0.35 


0.14 







ME, 


8.9 


6.8 


6.0 


5.2 


4.9 


4.7 


4.6 
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.0 2.0 



FIGURE 1. Volume differences, AV, versus normalized pressure, z, 
for the ME\ and other equations. 




1.0 2.0 



FIGURE 3. Volume differences, AV, versus normalized pressure, z, 
for the BE i and other equations. 




1.0 2.0 




1.0 2.0 



FlUURE 2. Volume differences, AV, versus normalized pressure, z, FIGURE 4. Volume differences, AV, versus normalized pressure, z, 
for the KE and other equations. for the 3DCE and other equations. 
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1.0 2.0 



Figure 5. Volume differences, AV, versus normalized /iressure, z, 
for the BE i and other equations. 

Although few actual experiments resulting in P-V 
values of appreciable accuracy extend past z ~ 0.5, the 
present exact, synthetic data curves are calculated up 
to p = 210 kbar, where z = 2.000. At this z value, 
V/'Vo = x~ i is of the order of 0.5 for these equations, 
being —0.64 for the 3DGE, for example. For 
z max = 2.000, A=58 p or z values, distributed roughly 
logarithmically, were used. For present purposes, 
larger z values were unnecessary. 

Clearly, AV curves for all pairs not involving the 
UTE, 3SE, or ME 2 will eventually reach a maximum, 
with AF max < 1, as z increases, then decrease toward 
zero since both Ps become arbitrarily small as z— » °°. 
As the figures show, the situation is different for the 
ME2 even within the present range. Since the param- 
eter values used here lead to V < for z > 1.85. AV 
values which involve ME2 volumes become arbitrarily 
large in magnitude as z increases beyond this point. 
Clearly, the ME 2 cannot be a useful model all the way 
to the point where it predicts zero or negative volumes. 
Nevertheless, it may be useful for a range ending 
sufficiently far below this point. 

Of what value are the results shown in figures 1 to 5? 
They are of considerable value because they show 
how well the various equations of state considered 
here may be discriminated under the best possible 
conditions. Suppose, for example, that we wish to 
discriminate between the KE and other equations and 
are able to measure volume only up to z = 0.1. Further, 
suppose that errors in p are negligible compared to 



those in V. Figure 2 then shows that to distinguish the 
3SE from the KE in the range 0<z<0.1, experi- 
mentally determined V values must be known to about 
one part in 10 4 , or to four decimal places, near z ~ 0.1. 
Even less uncertainty would be required for a smaller 
range. The BE 2 and KE cannot be reliably distin- 
guished without a precision of about three parts in 10" 
near 2: = 0.1 and higher precision for smaller z. Clearly, 
if the above precision has not been achieved, there 
would be no point in attempting to discriminate be- 
tween the equation pairs discussed for the data in 
question. Barsch and Chang [3] have discriminated 
between the BE 2 and KE for a situation where 
AV/Vo — 3 X 10 ~ 3 or more and have concluded that the 
BE^ was much better for their purposes than the KE. 
The present figure 2 results indicate that such dis- 
crimination is actually not significant with such pre- 
cision in AV, for the present set of parameter values, 
over a pressure range from zero up to at least 200 kbar. 

There are two reasons why we consider that the 
present curves represent the best possible discrimina- 
tion. First, there are always some random errors in the 
determination of pressure values. To first order, we 
may take the expected or "controlled" pressure values 
as exact and consider that the actual pressure errors 
are incorporated as additional random errors in the 
volume values. It is then this total volume error which 
must be used in determining whether the curves allow 
equation discrimination within ^ certain range of z. 
When parameter values are available, as from ultra- 
sonic measurements, they may be used in several 
equations of state to calculate exact volumes over a 
given z range. These volumes may then be directly 
compared with a set obtained by direct measurement. 
Clearly, if the total errors in the latter set are not 
appreciably smaller (over most or all of the z range) than 
the APs obtained with various equation pairs, no valid 
discrimination is possible. Even so, one of the several 
equations among which discrimination is impossible 
for the given z range may be far superior to the others 
for extrapolation beyond this range. Although all eight 
equations of figures 1 to 5 are indistinguishable for 
AV data of no better than 10~ 3 precision in the range 
0^z ^ 0.1, clearly there are important differences be- 
tween the predictions of the various equations for this 
same precision level at say z— 1.5. 

When an independently measured set of parameter 
values is unavailable, parameter value estimates must 
be obtained by fitting a model to the available data by 
some such procedure as least squares. Each different 
model fitted will then yield a different set of estimated 
parameter values. If AV values are obtained for a pair 
of models, using in each model the specific parameter 
values determined from a least squares fit of the data 
for the given model and range (case A), then the adjust- 
ment of the parameter values associated with the least 
squares procedure will generally lead to an appreciably 
different set of AV values than would have been ob- 
tained had the same parameter value set been used in 
each equation (case B). If the fits of the two equations 
for case A are sufficiently good, the corresponding AV 
values may nearly all be much smaller than those 
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obtained in case B with any single reasonable param- 
eter value set. But only one, at most, of the two sets of 
parameter values can be correct. Thus, one must 
proceed with extreme caution, and the small degree of 
discrimination possible from the case A fits and AVs 
is usually misleading. Further, any use of case A results 
outside of their fitting ranges is extremely dangerous. 

The most meaningful discrimination will be obtained 
from calculating AVs by the case B procedure, using 
the same most reasonable choice of parameter values 
in both equations. If the two equations under considera- 
tion seem to fit about equally well and no other param- 
eter value information is available, reasonable values 
to use in the case B discrimination are the averages of 
the two sets of values found from the least squares 
fittings. Because of the wide use of least squares 
procedures, these matters will be further discussed 
in the next section. 

The present case B results are closely related to 
some obtained by Barsch and Chang [3]. These authors 
compared, however, p-value predictions obtained from 
a certain lattice equation of state tailored for Csl 
with p values obtained from several other phenom- 
enological equations of state. They found, for example, 
that using the same parameter values the BE2 ap- 
proximated the lattice equation an order of magnitude 
better (in Ap) than did the KE. Although Barsch and 
Chang present V/Vo versus p curves for several of the 
equations of state considered herein, they do not give 
AV curves and are not primarily concerned with 
establishing what accuracy in V is needed, for a given 
p or z range to allow equation discrimination. 

Even though Barsch and Chang assert the superiority 
of the BE2 over the other phenomenological equations 
they considered, as already mentioned the BE 2 curve 
of the present figure 2 suggests that exceptionally ac- 
curate data or a very wide range will usually be re- 
quired to allow meaningful discrimination to be made 
between the BE 2 and the KE. Although Barsch and 
Chang's calculated | Ap| values for the BE2 and lattice 
equation were an order of magnitude smaller than those 
for the KE and lattice equation, the latter values them- 
selves were still considerably smaller for the range 
^ p ^ 200 kbar than either the errors in | Ap | calcu- 
lated from the BE 2 with experimental uncertainties in 
the parameters or those expected experimentally 
[3]. Thus, the actual discrimination between the BE 2 
and KE appears nugatory for this range. It seems 
doubtful that sufficiently accurate wide-range data yet 
exist to make adequate BE 2 -KE discrimination possible 
unless the situation is very different for appreciably 
different parameter values than those used here and 
those used by Barsch and Chang, an unlikely possibility. 

The curves of figures 1 to 5 are somewhat more 
general than they appear at first sight. First, since the 
normalized pressure variable z is used, the results are 
independent of the value of Kq. Second, since the Vo 
value used is quite close to unity, little specialization is 
introduced by the specific Vo value used. When V () 
differs appreciably from unity, the present curves may 
still be used with the AV values reinterpreted as 
AV/Vh values. For the UTE, ME,, and BE,, only the 



V 



ME 2 ' 



sign changes for the V BEl — Put 

V BEl — Vs SE curves shown in figure 5 indicate that the 
BEi remains a closer approximation to the other three 
equations over a wider range than if such changes of 
sign were absent. This result is perhaps one reason why 
the BEi has been found to be of relatively general 
applicability in the past. 

4. Least-Squares Comparisons 

Least squares procedures are frequently applied to 
noisy data for which the true underlying model is un- 
known and possibly nonlinear in some of the param- 
eters. Here we shall investigate the results of least 
squares fitting of exact data, especially with incorrect 
models. Such analyses, when the correct model and 
parameter values which generated the data are known, 



additional parameter 77 enters. This quantity is usually 
found to be within the range 3 < 77 < 8; thus, the 
present value, near 6, is fairly typical. Further, changes 
in 7) may be expected to change AF itself less than the 
Ps entering AV. On the other hand, the *P value used, 
near —7, is quite special since little is known thus far 
about the likely range of ^ for a variety of materials, 
and it probably can be positive as well as negative 
[1, 5]. Nevertheless, we suggest that the present curves 
may be used, at a zero to first order level, for an initial 
estimate of discrimination possibilities between various 
equations for other materials besides Csl at 150 °C. Of 
course, the next order of assessment would employ an 
estimated parameter set (Vo, K ih K' {) , and K[[ values) 
for the material in question. This set could then be 
used, as herein, to generate AV curves for comparison 
with the estimated total errors of the experiment, all 
incorporated into the V values. 

As examples of such zero-order assessment, we may 
consider the data of Monfort and Swenson [6], Kell and 
Whalley [7], and Vedam and Holton [8]. Monfort and 
Swenson studied potassium metal up to z ~ 0.4. Their 
volume data were given to four places, and they found 
a scatter of about 5 units in the last place. Although 
they primarily considered the BE,, the ME, was also 
introduced. The ME, curve of figure 5 shows a maxi- 
mum \AV\ for these equations of about 7X 10~ :i . When 
the Monfort-Swenson data is normalized to a Vo near 
unity, allowing comparison of V errors with present 
APs, one may estimate that the data are accurate to 
perhaps 3 X 10~ 3 in normalized volume. Comparison 
suggests that one might then just be able to distinguish 
between the BE, and ME, for this range and accuracy. 
One of the present authors [1] has considered dis- 
crimination between the 3DGE and 3SE for the °C 
water data of Kell and Whalley (z max ~ 0.05) and be- 
tween the 3DGE and ME 2 for the 50 °C water data of 
Vedam and Holton (z max ~ 0.44). Similar zero-order 
comparison of probable errors in V with the present 
AV curves suggests that the 3DGE-3SE discrimination 
was near the borderline of possibility and was probably 
not very meaningful, while the 3DGE~ME 2 discrimina- 
tion was somewhat more possible and certain. 

Finally, to the degree that the present AV curves are 
reasonably general, it is worth mentioning that the 
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can yield valuable information about the systematic 
errors arising from the use of the wrong model. Further, 
the use of exact data allows the usually mixed effects 
of random and systematic errors of this type to be 
entirely separated. Since figures 1 to 3 and 5 show that 
the 3DGE is, in some sense, close in its predictions to 
several of the other equations, it has been chosen here 
as the "correct" model for illustrative purposes. The 
exact data employed was thus generated by using in 
the 3DGE the 150 °C Csl parameter values already 
discussed. 

Table 3 shows the results of applying the least 
squares method in a few different situations of interest. 
Here and hereafter "linear" and "nonlinear" generally 
refer to linearity, or its absence, of the parameters 
entering the model. Thus, by a "linear" equation we 
will mean one linear in its parameters. The "linear" 
situation cited is actually rendered nonlinear in the 
parameters by the weighting used [lj. Even though the 
model is originally linear in the parameters, weighting 
of the independent variable will lead to nonlinear 
parameter dependence except in the special simple 
case (not considered here) of a linear relation between 
independent and dependent variables, hi a succeeding 
paper, we hope to investigate in some detail the im- 
portance of and degree of bias frequently arising in the 
A case of table 3 when random error is present. Here 
we continue to restrict attention to tin 1 exact data 
situation. 

The 3D(iE is written in table 1 in a form involving 
the parameters non linearly. This form was required by 
the constraint of using /£<>, 17, and M r as the basic param- 
eters in each equation listed in the table. On the other 
hand, the 3DGE may also be written in the linear form 



/>=£//,(*- iv 



(!) 



where A {) = when the Vq entering x ' — Vo/V is taken 
fixed and has its correct value (the procedure we 
shall use when Aw is a free parameter); A\ = Kol 
A 2 = (t?-1)£ () /2; and A 3 = (1/6) (tj 2 - 3t? + 2 + ¥)&. 
Clearly, direct linear least squares determination of the 
A; parameter estimates will allow corresponding V {) , 

Table 3. Possible errors in least squares estimates 





Parameter 


Parameter 


Situation 


estimates 


variance 
estimates 


A. Correct model, random 






errors in data 






Linear: correct weighting 


Unbiased, efficient 


Unbiased 


Linear: wrong weighting 


Unbiased, not 
efficient 


Biased 


"Linear": weighting of in- 


Biased 


Biased 


dependent variable 






Nonlinear 


Biased 


Biased 


B. Correct model, systematic 


Biased 


Biased 


errors in data 






C. Incorrect model 


Systematic 


Systematic 




error 


error 



/Co, tj, and ^ estimates to be calculated for comparison 
with their true values. Further, comparison of cor- 
responding nonlinear and linear least square parameter 
estimates will allow bias arising from nonlinear least 
squares to be indentified and quantified. 

The following definitions are useful in comparing 
least-squares parameter values with exact values. Let 
be a specific parameter of the model: then denote the 
true value of 6 (here known) by 0o and the least-squares 
estimate as 0. The relative 1 error of the estimate is then 
8 = (0 — 6{)) I do. When no random errors are present. 
8i will measure the systematic error in the / 1 1 1 param- 
eter value, It is also of interest to compare the param- 
eter error (0 — O {) ) with the standard deviation (s ( i)h 
obtained for a given least squares estimate of 6. To do 
so, we define A = | (5- 0„)/(s,/)« | = |ft,8/(s rf )« |. This 
measure will indicate possible systematic error in 

We have been discussing least squares results in the 
above as though they were exact solutions of the leasl 
squares equations. It is not widely appreciated that 
all the usual least squares computer routines may yield 
very inaccurate parameter values because of round-off 
errors | ( )|. For example, if Gaussian elimination with 
pivoting is used to solve the least squares equations, 
the number of accurate decimal digits in a ft, A, is 
A ~ (C — n+ 1± 1), where C is the number of (equiva- 
lent) decimal digits carried in the computer calculation 
and // is the number of free parameters. Clearly, if 
m>C, results of little value are likely to be obtained. 
Expression for A of this type were originally derived 
for linear least squares fitting of polynomials, but they 
seem to apply at least approximately to nonlinear equa- 
tions as well. Recently, \\ ampler [ 1()| has made a more 
detailed study of the matter for polynomials and dis- 
cussed more complex routines which can yield very 
substantially higher solution accuracy. 

The effects of roundoff are illustrated by the results 
of table 4. Elimination with pivoting was used to carry 
out least squares fitting of the 14-figure 3DGE data using 
the 3DGE equation in both its linear and nonlinear 
forms. Since c=14 and ti=4, ^~11±1. In Table 4, 
the 5/ are calculated using rV =K , 0i = /Co, 0- 2 = rj, 
and 0.) = ^. The quantity Sd is the standard deviation 
for the fit itself. The results show values of A between 
about 14 and slightly less than 11, in rough agreement 

TABLE 4. Least squares results in the absence of 
systematic error: exact 3DGE data fitted by the 
3DGE model 
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Linear < 


equation 


Nonlineai 


• equation 




p-Weighting 


^-Weighting 


p-Weighting 


^-Weighting 


So 

8, 
8., 
8., 


-2.7 x 10- l: <) 

3.3 x 10- 13 

-1.5 X10- 12 

-1.9x10-" 


04o = 

-2.4X 10 1:{ ) 
3.0 x 10- 13 
-1.3X 10 '-' 

-1.7X 10-" 


-7. Ox 10- 15 

3.5 x 10 1:{ 
-1.4X 10 12 
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-7.0X 10" 

2.9 x lO" 13 
-1.1 X 10- 12 
-1.5X 10-" 


Sd 


2.3 x 10-' 2 


2.0X 10 14 


2.3X 10-' 2 


2.0 x 10-" 
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with the formula. There appears to be no significant 
tendency for the linear results to be better than the 
nonlinear ones, and one can scarcely conclude that 
much of the bias of table 3 is showing up here. In fact, 
bias is only important when random errors are present; 
in the A cases of table 3, bias approaches zero as the 
random error goes to zero. Incidentally, since 80 is 
zero in the linear case when V {) is taken fixed at its 
exact value, A {) is given in its place; since the true 
value of A{) is zero, this is an absolute, not relative error. 

The results of table 4 were calculated with A^ = 37 
points, covering the range ^ z ^ 0.476. Let the maxi- 
mum value of z be denoted z r . In earlier work [1], 
weighting of both the dependent and independent 
variable data values was discussed. The related stand- 
ard deviations were denoted en for the V variable and 
(T/j for the p variable. The p-weighting of table 4 takes 
(T p = l and cr\=0 (weighting of dependent variable 
only), while the V- weighting uses o> = 0, oy = 1. In the 
linear equation case, the F-weighting chosen leads to 
somewhat different weighting of the actual independent 
variable x used [1 ]. Table 4 indicates slightly improved 
results for the F-weighting over the p-weighting, and 
no bias arising from V- weighting in the nonlinear situa- 
tion is apparent. The differences between the s/s for 
p and V weighting arise because the p-weighting sa is 
a measure of the least-square fit residuals when they 
are all in pressure and is here in kilobars, while for 
F-weighting the residuals are all forced to be in re- 
duced volume and s ( \ is then dimensionless. The ratio 
between the s/s is roughly Ko. 

Although we shall use the usual inaccurate Gaussian- 
elimination-with-pivoting solution of the least squares 
equations in the following, all inaccuracies introduced 
thereby are four or more orders of magnitude smaller 
than the systematic errors we consider. Such syste- 
matic errors in parameters and standard deviations are 
illustrated in table 5, where the 3DGE data are fitted 
with p-weighting using the incorrect BE 2 model. Re- 
sults for 8, and A; are first given for two different ranges 
of 2, from zero to ~ 0.048 and ~ 0.48. Note the strong 
increases in these error measures both with range and 
with the index i. Also included in the table are fitting 
results obtained for the complement range, all points 
contained in the second range but not in the first. As 
might be expected, the parameter estimates are 
somewhat worse for this coverage than for the largest 
span shown, even though s ( i itself is somewhat better. 



Table 5. Least squares results showing systematic 
errors: exact 3DGE data fitted by the BE 2 model 
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3.4x10' 


9.7 
14.0 
24.0 
45.0 


Sd 


1.31 X10- 8 


1.02X10 - 3 


4.62 XIO" 4 



All nonlinear least squares fitting in the present 
paper has been carried out using the Deming itera- 
tive method of solution [1, 11 J. Although this is an 
approximate method [12], the resulting errors in the 
estimated parameter values are generally negligible. 
O'Neill et al. [12] have presented a more accurate 
iterative solution of the least squares problem with 
weighting of both dependent and independent variable, 
applicable only for polynomial equations. We have 
recently generalized and improved this solution so 
that it applies to equations of any form and converges 
much more rapidly [13]. This method, applied to the 
situations of table 5, yields essentially the same 6\- 
values as those in the table but A,-'s some 25 to 40 
percent larger than the tabular values. These increases 
thus mainly arise from smaller (s<fVs produced by the 
new least squares solution. Although the new method 
leads to an essentially exact (in the sense of iterative 
convergence) least squares solution when round-off 
errors are negligible, the results cited above and those 
in the table show the presence of large systematic 
errors in 8, and A, arising from wrong model choice. 
The rest of the present paper is primarily concerned 
with oYs and p-weighting, where the differences 
between the Deming and improved estimates are 
negligible. 

Although table 5 gives one some idea of syste- 
matic error effects, much more is provided by the least- 
squares results of figures 6 to 12. The same exact 
3DGE data were fitted with the various other equations 
for four different ranges, all including z = 0. The 
four values of z r used were -0.048, 0.143, 0.476, 
and 2.00, and the corresponding number of z values 
were, respectively, 19, 29, 37, and 58. All points used 
in a given fitting were included in the ones with larger 

Z r . 
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Figure 6. Parameter relative errors, 8 i9 and standard deviation of 
fit, s d , versus fitting range (0 «S z «S z r ) for least squares fitting of 
exact 3DGE data with the BE\. 
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FIGURE 7. Parameter relative errors, h x and standard deviation of 
Jit, s d , versus fitting range (()=£ z ^ z T ) for least squares fitting of 
exact 3DGE data with the ME X . 
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Figure 9. Parameter relative errors, ftj, and standard deviation of 
fit, S,/, versus fitting range (0 ^ z =S z T ) for least squares fitting of 
exact 3DGE data with the BE 2 . 
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Figure 8. Parameter relative errors, 8 X , and standard deviation of FIGURE 10. Parameter relative errors, 8j, and standard deviation of 
fit, Sd, versus fitting range (0 =S z ^ z r ) for least squares fitting of fit, s d , versus fitting range (0 =^ z ^ z r )/br least squares fitting of 
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FIGURE 11. Parameter relative errors, 8,, and standard deviation of 
fit, s d , versus fitting range (0 ^ z =£ z r ) for least squares fitting of 
exact 3DGE data with the KE. 

Figures 6 to 12 show how the systematic errors, 
represented by the oYs and by s,/, change as the fitting 
range is extended. As usual, dashed lines indicate 
negative values. All 8/ results shown were obtained 
with p-weighting: crr = 0, <t p =1. Results obtained 
with F-weighting were closely similar. Unlike the 8,'s, 
which are relative, the s/s are absolute and, with 
p-weighting, measure the overall goodness of fit in 
pressure units, as already mentioned. Thus, for 
example, figure 6 indicates that Sd for the BE] fitting 
over the range 0^z^2 is nearly 0.1 kbar. All Sd 
curves were obtained with p-weighting except the one 
marked (Sd)v on figure 9. Here we compare the s ( fs 
obtained from p and V weighting. The (sd)v values are 
somewhat more than K {) times smaller than (s ( t) 1} 
values here. Note that, as expected, (s (1 )v and 6 () , 
the relative error in Viu are quite close together over 
much of the range. 

For V weighting. Sd is an overall measure of the 
residuals in V. Its absolute value in figure 9 at z = 0.143 
of about 3 X 10" 7 (the maximum magnitude of a volume 
residual was ~ 6.5 X IO -7 ) is about two orders of magni- 
tude smaller than the AF '= V BE2 ~V DGE of ~5X IO" 5 
shown in figure 3 for the same z. But this last figure is 
that applying when the correct parameter values are 
used in both equations. As expected, the least squares 
parameter adjustment in the BEo fitting of the exact 
3DGE data makes it difficult to conclude (without in- 
dependent knowledge of parameter values) that the 
BE 2 is the wrong model, as it is here. With some ran- 
dom errors in the 3DGE model data, least squares 
fitting using the BE 2 and KE, for example, would again 
generally lead to results which wouldn't allow one to 



identify either the BE 2 or KE as an incorrect model, 
even though they both would be. 

The results of figures 6 to 12 show that when the 
range is extended, relative errors in all the parameters 
increase when wrong models are used. Further, the 
higher-order parameters are more inaccurate than the 
lower-order ones for all the ranges shown. Not much 
added accuracy in the higher-order parameters can 
be obtained by reducing the range and, in practical 
cases where random error is present, generally no 
added accuracy will be achieved by such reduction. 

Figure 10 stops with a z r of 0.476 because the volume 
predicted by the ME 2 is negative for z ^ 1.85, preclud- 
ing a meaningful fit with z r — 2.00. Note that 8 ;? for the 
KE and 3SE is so large that its values must be divided 
by 10 and 100, respectively, to allow plotting on the 
present scale. For the 3SE, even 8 2 must be divided by 
10 as well. These results illustrate an important gen- 
eral point. The figures show that the BE 2 and KE are 
the best least-squares simulators of the 3DGE model 
as far as s ( i is concerned. Yet even for the relatively 
low z value of 0.143 (p= lS^kbar), |^| is about 10 per- 
cent high for the BE 2 and "9 is of even the wrong sign 
for the KE. The average residuals arising from syste- 
matic error would, when all in volume, be mostly less 
than IO 6 in magnitude here. Even for the best data 
currently available such small residuals would be ob- 
scured by random error. Thus we see that it is possible 
that two different equations, both wrong (as here) or 
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FIGURE 12. Parameter relative errors, 8 i? and standard deviation of 
fit, s d , versus fitting range (0 =S z ^ z r ) for least squares fitting of 
exact 3DGE data with the 3SE. 
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one wrong and one correct, may not be distinguishable 
by goodness of fit criteria, yet one may predict far 
better parameter values than the other. In the absence 
of other information, such as firm knowledge of the 
correct model or independent determinations of some 
of the parameter values, one will evidently always 
stand an appreciable chance of picking a "best" model 
which yields some quite poor parameter estimates. 
The better the accuracy of the data and the wider its 
range, the better the higher-order parameter estimates 
will be since the final model chosen will be forced to 
be closer to the correct model to achieve an adequate 
fit. 

The monotonic increase of 8, and s ( i with fitting range 
illustrated in figures 6 to 12 is, of course, indicative of 
the use of an incorrect and inadequate model and is 
by no means limited to the equation of state area. In 
most if not all cases of practical interest, we may expect 
to find this sort of behavior: the wider the range used 
in least-squares fitting of a possible, "close," but still 
incorrect model, the greater will be s ( i and the param- 
eter error magnitudes. It should, however, be remarked 
that this conclusion only applies in the usual case 
where the model is not asymptotically correct as the 
range is extended indefinitely. The wider the range 
used, generally the more difficult it will be for an 
incorrect model to simulate the correct one. 

This increase of errors with range may frequently be 
used in practical cases as a powerful means of dis- 
criminating against incorrect models. When random 
errors in the data are sufficiently small that the sys- 
tematic errors arising from an incorrect model choice 
dominate s r /, it will generally increase with the fitting 
range, as illustrated here. Such an increase thus clearly 
signals an incorrect model choice for the range of data 
fitted. Since most models only apply adequately in any 
case over a limited range of a variable such as pressure 
or temperature, extension of the fitting range beyond 
the region of applicability of the best available model 
will always eventually result in an increase in Sd> Thus, 
in any least squares fitting where the range of applica- 
bility of the model used in unknown, extrapolation 
outside the fitted range of data should be approached 
with the utmost caution and avoided if possible. 

The present paper deals with exact data and actual 
relative errors of parameters, but true parameter errors 
will not be available in a usual experimental situation. 
Nevertheless, when s<t increases because of the wrong 
model choice, the estimated parameter standard de- 
viations will generally increase for the same reason. 
Thus, these quantities, ordinary results of a least 
squares fitting, may also be used along with s f f to help 
discriminate against an inadequate model. 

There are some interesting general aspects to the 
present results obtained with least-squares fitting of 
the wrong model. The residuals (here given by observed 
values minus calculated values) show the following 
behavior: The number of runs (number of sign changes 
plus one), u, for the MEi, BEi, and UTE, for which 
n = 3, is 4, while u= 5 for the remaining equations for 
which n = 4. The general result, u — (n + 1 ) , is not very 
surprising but bears emphasizing. Further, the sign of 
the first residual run (which, together with knowledge 



of u determines the signs and order of all the runs) is 
specific to the equation considered. For the present 
fitting of 3DGE data, this sign is -f , — , — , + , — , — , + for 
the ME!, BE,, UTE, ME 2 , KE, BE>. and 3SE, respec- 
tively. The number of runs and their sign distributions 
were invariant in the present situation to the following: 
(a) p or V weighting, (b) the range of the data and its 
placement (all low p, all high />, all in the middle, etc.), 
and (c) the sign of ^V . Even though not all extremes were 
investigated, this high degree of pattern invariance is 
likely to be quite general and may itself be of con- 
siderable usefulness in helping to distinguish models. 

Although we have not done it, one could readily 
establish a matrix of first signs obtained using data 
calculated from one of the present eight specific equa- 
tions and fitted with another one of the eight. Then, 
in practical situations where it was believed that the 
correct model was one of the eight, many possibilities 
could be quickly eliminated by comparison with the 
sign of the first residual run obtained on fitting the 
actual data. This would only work, of course, provided 
random errors were considerably smaller than sys- 
tematic ones and hence didn't appreciably perturb the 
residual pattern. With many data points, considerable 
perturbation of this kind could be tolerated, however, 
since decisions could be made on the basis of a 
smoothed residual pattern rather than tin 1 actual noisy 
pattern. 

A partial comparison of the above type has been 
made earlier for the ME, and UTE [14]. There, V {) 
was taken fixed, so n = 2. As expected, u was found to 
be three for both UTE fitting of exact ME, data and for 
ME, fitting of UTE data. The initial run signs were 
+ , — , respectively, for the above two fittings. 

5. Summary 

This paper has been primarily concerned with dis- 
cussing methods of discriminating between specific 
equations of state and has demonstrated considerable 
limitations on the possibility of adequate discrimina- 
tion between "close" equations. We have found the 
somewhat surprising result that equations which can- 
not be adequately discriminated on the basis of least 
squares goodness of fit over even a wide pressure range 
may yet lead to estimates of higher-order parameter 
relative errors differing in sign and by an order of 
magnitude in absolute value for even a narrow pressure 
range, much less a wide one. The present methods, 
results, and conclusions can be generalized to a con- 
siderable degree to apply to model discrimination out- 
side the equation of state area and are pertinent for 
linear models and for those nonlinear in their param- 
eters, independent variable, or both. Thus, the follow- 
ing general conclusions, based on the present specific 
results, are likely to apply widely to the general data 
analysis field. 

More than one mathematical model should usually 
be tested against the data in order to select, if possible, 
that model which fits best by objective criteria. As the 
range of data is progressively increased for which least 
squares fittings are carried out, the initial or eventual 
appearance of increases in s ( i and (sd)u, indicates the 
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FIGURE 13. General block diagram for data analysis. 

presence of systematic fitting error arising from an 
inadequate model choice. Such error will also usually 
show up in highly correlated residuals exhibiting, at 
least approximately, a number of sign-changes equal 
to the number of fitted parameters. The range of a 
causal experimental variable such as pressure, voltage, 
temperature, etc. should be increased to the maximum 
degree possible in order to allow the testing of a model 
for adequacy over the widest practical data range. 

When two or more models have been found that 
represent the data over a given range with approxi- 
mately the same goodness of fit and without signs of 
systematic errors from wrong model choice, it is still 
possible that one or more models may yield much 
better or much worse least-squares parameter esti- 
mates than the others. Additional independent infor- 
mation about likely parameter value ranges will usually 
then be necessary to allow a selection of the most 
appropriate equation for parameter estimation. Extrap- 
olation of a given model-parameter value set beyond 
the range of data on which it is based is always 
dangerous. 

When data smoothing or interpolation is the object, 
the possibility of discrimination between two models 
which yield equally good least squares fits to the data 
should be examined by the case B procedure of section 
3. If the differences in dependent variables calculated 
with the same reasonable set of parameter values in 
each model are comparable to or smaller than the 
estimated random errors in the data, discrimination is 
impractical for that data set. 

Figure 13 shows, in very diagrammatic form, ap- 
propriate steps in data analysis aimed at establishing 
a "best" model (including specific parameter value 
estimates). Some of the actually interrelated steps 
involved in this figure are presented differently in the 
flow chart of figure 14. This diagram is included for the 
benefit of those readers who may wish to apply the 
procedures discussed in this paper to other discrimi- 
nation and parameter estimation problems. 

For figure 14 we have assumed that a data set over 
a range, R max ^ nas been taken, and that we wish to test 
potential models over the maximum range if possible. 



The flow chart orders the tests as (1) case B, (2) runs, 
and (3) case A. If no models appear appropriate after 
the first series of tests, provision is made for decreas- 
ing the range of the data used in the tests in order to 
determine the acceptable maximum range for param- 
eter estimation. 

In the flow chart, we have abbreviated the test for 
case B by the notation |Ay/j | <o> Here we mean that 
all or nearly all of the deviations should be less than 
the expected errors in the data. Note that "nearly all" 
is appropriate because of the possible presence of 
random outliers. For the same reason, the test u>n+2 
should also be considered approximate and applied 
judiciously. Note also that <r y may vary with x (het- 
eroscedastic case); the test should be so applied when 
appropriate. Other symbols introduced here are e s , 
defined to be the acceptable level for standard devia- 
tion of the least squares fitting, and e f /, defined as the 
level below which standard deviations of two separate 
fits are indistinguishable. 

Good data are usually expensive, yet too little 
adequate data analysis is the general rule. It is better 
to do too much such analysis than too little. 
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Figure 14. 



Flow diagram for discrimination and parameter estima- 
tion tests. 
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